﻿#ifndef SWEEP_CIRCLE_VORONOI
#define SWEEP_CIRCLE_VORONOI

struct Event;
struct Ellipse;
struct EllipseArc;
struct ActiveEdge;
struct EllipseArcNoPt;
struct ActiveEdgeNoPt;



template<class T>
class __align__(16) Pool{
public:

	__device__ void initial(T* origin,int size,int capacity){
		d_Pool = origin;
		_size = size;
		_capacity = capacity;
	}
	__device__ bool insert(const T& tmpEle){
		if( _size >= _capacity ){
			cuPrintf("insert pool full!!\n");
			return false;
		}
		d_Pool[_size++] = tmpEle;
		return true;
	}
	__device__ bool alloc( T* & tmpPt){
		if( _size >= _capacity){
			cuPrintf("alloc pool full!!!\n");
			return false;
		}
		tmpPt = &d_Pool[_size++];
		return true;
	}
	__device__ int size(){
		return _size;
	}
	__device__ T& operator[](const int& pos){
		return d_Pool[pos];
	}
//private:
public:
	T* d_Pool;
	int _size;
	int _capacity;
};


/* 32
	ellipse 24
*/

struct __align__(32) Ellipse
{
	float2 majorFocus;
	float2 minorFocus;
	float rotateAngle;
	float c;
	__device__ void initial(float2 _majorFocus,float2 _minorFocus){
		majorFocus = _majorFocus;
		minorFocus = _minorFocus;
		rotateAngle = atan2f(minorFocus.y - majorFocus.y,minorFocus.x - majorFocus.x);
		c=sqrtf((majorFocus.x-minorFocus.x)*(majorFocus.x-minorFocus.x)+
			(majorFocus.y-minorFocus.y)*(majorFocus.y-minorFocus.y))/2;

	}

	__device__ __host__ float GetA(float d){
		return d / 2;
	}
	__device__ __host__ float GetC(){
		return c;
	}
	__device__ __host__ float GetE(float d){
		return GetC()/GetA(d);
	}
	__device__ __host__ float GetP(float d){
		float a = GetA(d);
		float c = GetC();
		return a * a / c - c;
	}
	__device__ __host__ float GetRotateAngle(){
		return rotateAngle;
	}
	__device__ __host__ float2 GetPointAtAngleParam(float d,float angle){
		float e = GetE(d);
		float p = GetP(d);
		if( !(p>WXN_EPSILON) ){
			//assert(p>WXN_EPSILON);
			//cuPrintf("p>WXN_EPSILON ,error! c %f e %f d %f p %f\n" , c , e , d , p);
			
			return make_float2(FLT_MAX,FLT_MAX) ;
		}
		float alpha = GetRotateAngle();
		float rou = e * p / ( 1 - e * cosf(angle - alpha) );
		float2 pt;
		pt = make_float2(rou * cosf(angle),rou * sinf(angle));
		pt.x += majorFocus.x;
		pt.y += majorFocus.y;
		return pt;
	}
	__device__ __host__ static float2 GetIntersectionAngleParam(Ellipse& ellipse1,Ellipse& ellipse2,float d){
		float e1 = ellipse1.GetE(d);
		float e2 = ellipse2.GetE(d);
		float p1 = ellipse1.GetP(d);
		float p2 = ellipse2.GetP(d);
		float alpha1 = ellipse1.GetRotateAngle();
		float alpha2 = ellipse2.GetRotateAngle();
		float a = p2 * cos(alpha1) - p1 * cos(alpha2);
		float b = p2 * sin(alpha1) - p1 * sin(alpha2);
		float c = (e2 * p2 - e1 * p1) / (e1 * e2);
		float len = sqrtf(a*a+b*b);
		a /= len;
		b /= len;
		c /= len;
		float t = atan2f(b,a);
		float t2;
		if( c >= 1 )
			t2 = 0;
		else if ( c <= -1)
			t2 = M_PI;
		else 
			t2 = acosf(c);
		return make_float2(t-t2,t+t2);
	}
};



struct __align__(64) EllipseArc : Ellipse
{
	
	bool outofdate;
	ActiveEdge* leftActiveEdge;
	ActiveEdge* rightActiveEdge;

	__device__ void initial(float2 _majorFocus,float2 _minorFocus,ActiveEdge* _leftActiveEdge,ActiveEdge* _rightActiveEdge){
		Ellipse::initial(_majorFocus,_minorFocus);
		leftActiveEdge=_leftActiveEdge;
		rightActiveEdge = _rightActiveEdge;
		outofdate = false;
	}

	__device__ void GetVanishEvent(Event& evt,Pool<float2>& vertexPool );
	std::vector<std::vector<float2> > EllipseArc::GetSegments(float d,float2 tmpCenter,float2 tmpMinPoint,float2 tmpMaxPoint);

	__device__ bool Zero(float d){
		float2 angles = GetAngles(d);
		return angles.y - angles.x < WXN_EPSILON;
	}
	__device__ float2 GetAngles(float d);


};



struct __align__(32) ActiveEdge{
	EllipseArc* leftArc;
	EllipseArc* rightArc;
	int fixedVertex;
	float2 direction;
	__device__ bool initial(EllipseArc* _leftArc,EllipseArc* _rightArc,int _fixedVertex){
		leftArc = _leftArc;
		rightArc = _rightArc;
		fixedVertex = _fixedVertex;
		direction = GetDirection(_leftArc->minorFocus,_rightArc->minorFocus);
		
		if( direction.x == FLT_MAX && direction.y == FLT_MAX) return false;

		return true;
	}
	/*ActiveEdge(EllipseArc* _leftArc,EllipseArc* _rightArc,int _fixedVertex):leftArc(_leftArc),rightArc(_rightArc),fixedVertex(_fixedVertex){
		float2 p1 = _leftArc->minorFocus;
		float2 p2 = _rightArc->minorFocus;
		direction = GetDirection(p1,p2);
	}*/
	__device__ __host__ float2 GetCurrentPos(float d ){
		if( leftArc == NULL || rightArc == NULL ){
			//cuPrintf("error in GetCurrentPos\n");
			return make_float2(0,0);
		}
		if( leftArc->GetP(d) < WXN_EPSILON && rightArc->GetP(d) < WXN_EPSILON ){
			return leftArc->majorFocus;
		}
		if(leftArc->GetP(d) < WXN_EPSILON){
			float angle = atan2f(leftArc->minorFocus.y - leftArc->majorFocus.y,leftArc->minorFocus.x - leftArc->majorFocus.x);
			return rightArc->GetPointAtAngleParam(d,angle);
		}
		if(rightArc->GetP(d) < WXN_EPSILON){
			float angle = atan2f(rightArc->minorFocus.y - rightArc->majorFocus.y,rightArc->minorFocus.x - rightArc->majorFocus.x);
			return leftArc->GetPointAtAngleParam(d,angle);
		}
		float2 params = Ellipse::GetIntersectionAngleParam(*leftArc,*rightArc,d);
		float2 p1 = leftArc->GetPointAtAngleParam(d,params.x);
		float2 p2 = leftArc->GetPointAtAngleParam(d,params.y);
		if((p2.x-p1.x)*direction.x + (p2.y-p1.y)*direction.y>0){
			return p2;
		}else{
			return p1;
		}
	}
	float2 __device__ GetDirection(float2 p1 , float2 p2)
	{
		if( p1.x == p2.x && p1.y == p2.y ){
			cuPrintf("direction error!\n");
			return make_float2(FLT_MAX,FLT_MAX);
		}
		float x = p1.y - p2.y;
		float y = p2.x - p1.x;
		float len = sqrt(x*x+y*y);
		return make_float2(x/len,y/len);
	}
};
struct __align__(32) ActiveEdgeNoPt{
	int leftArcPos;
	int rightArcPos;
	int fixedVertex;
	float2 direction;
	__device__ void initial(ActiveEdge* activeEdge,EllipseArc* arcStartPos){
		fixedVertex = activeEdge->fixedVertex;
		direction = activeEdge->direction;
		if( activeEdge->leftArc == NULL ){
			leftArcPos = -1;
		}else{
			leftArcPos = activeEdge->leftArc - arcStartPos;
		}
		if( activeEdge->rightArc == NULL ){
			rightArcPos = -1;
		}else{
			rightArcPos = activeEdge->rightArc - arcStartPos;
		}
	}
};

struct __align__(64) EllipseArcNoPt:Ellipse{
	bool outofdate;
	int leftActiveEdgePos;
	int rightActiveEdgePos;
	__device__ void initial(EllipseArc* arc,ActiveEdge* edgeStartPos){
		majorFocus = arc->majorFocus;
		minorFocus = arc->minorFocus;
		rotateAngle = arc->rotateAngle;
		c = arc->c;
		outofdate = arc->outofdate;
		if( arc->leftActiveEdge == NULL ){
			leftActiveEdgePos = -1;
		}else{
			leftActiveEdgePos = arc->leftActiveEdge - edgeStartPos;
		}
		if( arc->rightActiveEdge == NULL ){
			rightActiveEdgePos = -1;
		}else{
			rightActiveEdgePos = arc->rightActiveEdge - edgeStartPos;
		}
	//	cuPrintf("leftactiveEdgepos %d rightActiveedgepos %d\n" , leftActiveEdgePos , rightActiveEdgePos );
	}
}; 



enum EventType {NewEllipse,EllipseVanish};

struct __align__(16) Event{
	float distance;
	int type;
	float2 site;
	EllipseArc* arc;
	ActiveEdge* leftActiveEdge;
	ActiveEdge* rightActiveEdge;
	__device__ __host__  bool operator<(const Event& other)const{
		return distance < other.distance;
	}
};
struct __align__(16) EventNoPt{
	int arcPos;
	int leftActiveEdgePos;
	int rightActiveEdgePos;
	__device__ void initial(Event* evt,EllipseArc* arcStartPos,ActiveEdge* activeEdgeStartPos ){
		if( evt->arc == NULL ){
			arcPos = -1;
		}else{
			arcPos = evt->arc - arcStartPos;
		}
		if( evt->leftActiveEdge == NULL){
			leftActiveEdgePos = -1;
		}else{
			leftActiveEdgePos = evt->leftActiveEdge - activeEdgeStartPos;
		}
		if( evt->rightActiveEdge == NULL ){
			rightActiveEdgePos = -1;
		}else{
			rightActiveEdgePos = evt->rightActiveEdge - activeEdgeStartPos;
		}

	}
};


__device__ void EllipseArc::GetVanishEvent(Event& evt,Pool<float2>& vertexPool ){
	evt.type = EllipseVanish;
	
	evt.distance = FLT_MAX;
	evt.arc = this;
	if(leftActiveEdge==NULL || rightActiveEdge == NULL){
		return;
	}
	if(leftActiveEdge->fixedVertex == rightActiveEdge->fixedVertex){
		return;
	}
	if(fabsf(leftActiveEdge->direction.x * rightActiveEdge->direction.x+leftActiveEdge->direction.y * rightActiveEdge->direction.y)>1.0 - FLT_EPSILON){
		return;
	}
	float e1 = vertexPool[leftActiveEdge->fixedVertex].x - vertexPool[rightActiveEdge->fixedVertex].x;
	float e2 = vertexPool[leftActiveEdge->fixedVertex].y - vertexPool[rightActiveEdge->fixedVertex].y;
	float a = -leftActiveEdge->direction.x;
	float b = rightActiveEdge->direction.x;
	float c = -leftActiveEdge->direction.y;
	float d = rightActiveEdge->direction.y;
	float mod = a * d - b * c;
	float tLeft = d * e1 - b * e2;
	tLeft /= mod;
	float tRight = -c * e1 + a * e2;
	tRight /= mod;
	if(tLeft < -FLT_EPSILON || tRight < -FLT_EPSILON){
		return;
	}
	if(tLeft<0)
		tLeft = 0;
	if(tRight<0)
		tRight = 0;
	float2 circumcenter;
	circumcenter.x = vertexPool[leftActiveEdge->fixedVertex].x + tLeft * leftActiveEdge->direction.x;
	circumcenter.y = vertexPool[leftActiveEdge->fixedVertex].y + tLeft * leftActiveEdge->direction.y;
	float dis = sqrtf((circumcenter.x - minorFocus.x)*(circumcenter.x - minorFocus.x)+(circumcenter.y - minorFocus.y)*(circumcenter.y-minorFocus.y));
	dis += sqrtf((circumcenter.x - majorFocus.x) * (circumcenter.x - majorFocus.x)+(circumcenter.y-majorFocus.y)*(circumcenter.y - majorFocus.y));
	evt.distance = dis;
	evt.leftActiveEdge = leftActiveEdge;
	evt.rightActiveEdge = rightActiveEdge;
	evt.site = circumcenter;
	return;
}

__device__ __host__ float2  EllipseArc::GetAngles(float d)
{
	if(leftActiveEdge==NULL){

		//cuPrintf("getangles 1\n");
		return make_float2(0,2*M_PI);
	}
	if(rightActiveEdge == NULL ){
		//cuPrintf("error in GetAngles");
		return make_float2(0,0);
	}
	if(GetP(d)<FLT_EPSILON){
		if(fabsf(minorFocus.y - majorFocus.y) + fabsf(minorFocus.x - majorFocus.x)<WXN_EPSILON){
			//cuPrintf("error in get angles\n");
			return make_float2(0,0);
		}else{
			//cuPrintf("getangles 2\n");
			return make_float2(GetRotateAngle(),GetRotateAngle());
		}
	}
	float2 leftP = leftActiveEdge->GetCurrentPos(d);
	float2 rightP = rightActiveEdge->GetCurrentPos(d);
	if( fabsf(rightP.y - majorFocus.y)+fabsf(rightP.x - majorFocus.x)< WXN_EPSILON ){
		//cuPrintf("error in get angles\n");
		return make_float2(0,0);
	}
	float angle1 = atan2f(rightP.y - majorFocus.y ,rightP.x - majorFocus.x);
	if( fabsf(leftP.y - majorFocus.y) + fabsf(leftP.x - majorFocus.x ) < WXN_EPSILON ){
		//cuPrintf("error in get angles\n");
		return make_float2(0,0);
	}
	float angle2 = atan2f(leftP.y - majorFocus.y,leftP.x - majorFocus.x );
	if( angle2 < angle1 - WXN_EPSILON ){
		angle2 += 2 * M_PI;
	}else if ( angle2 < angle1 + WXN_EPSILON && leftActiveEdge->leftArc == rightActiveEdge->rightArc ){
		angle2 += 2 * M_PI;
	}else if( angle2 - angle1 > 2 * M_PI - WXN_EPSILON){
		angle1 = angle2;
	}
	//cuPrintf("getangles 3\n");
	return make_float2(angle1,angle2);
}

std::vector<std::vector<float2> > EllipseArc::GetSegments(float d,float2 tmpCenter,float2 tmpMinPoint,float2 tmpMaxPoint)
{
	float minAngle = -M_PI;
	float maxAngle = M_PI;
	int xRegionCode = 0;
	int yRegionCode = 0;
	if ( tmpCenter.x > tmpMaxPoint.x )
	{
		xRegionCode = 1;
	}
	else if ( tmpCenter.x < tmpMinPoint.x )
	{
		xRegionCode = -1;
	}
	if ( tmpCenter.y > tmpMaxPoint.y )
	{
		yRegionCode = 1;
	}
	else if (tmpCenter.y < tmpMinPoint.y)
	{
		yRegionCode = -1;
	}

	if (xRegionCode == 0 && yRegionCode == 0)
	{
		minAngle = -M_PI;
		maxAngle = M_PI;
	}
	else if (xRegionCode == 1 && yRegionCode == 0)
	{
		minAngle = atan2(tmpMaxPoint.y - tmpCenter.y, tmpMaxPoint.x - tmpCenter.x);
		maxAngle = atan2(tmpMinPoint.y - tmpCenter.y, tmpMaxPoint.x - tmpCenter.x) + 2 * M_PI;
	}
	else
	{
		float angles[4] = {atan2(tmpMaxPoint.y - tmpCenter.y, tmpMaxPoint.x - tmpCenter.x),
			atan2(tmpMaxPoint.y - tmpCenter.y, tmpMinPoint.x - tmpCenter.x),
			atan2(tmpMinPoint.y - tmpCenter.y, tmpMinPoint.x - tmpCenter.x),
			atan2(tmpMinPoint.y - tmpCenter.y, tmpMaxPoint.x - tmpCenter.x)};
		minAngle = 10;
		maxAngle = -10;
		for (int i = 0; i < 4; ++i)
		{
			if (angles[i] > maxAngle)
				maxAngle = angles[i];
			if (angles[i] < minAngle)
				minAngle = angles[i];
		}
	}
	
	std::vector<std::vector<float2> > resPoints;
	if (GetP(d) < FLT_EPSILON)
	{
		resPoints.push_back(std::vector<float2>());
		resPoints[0].push_back(minorFocus);
		if (leftActiveEdge == NULL)
		{
			resPoints[0].push_back(majorFocus);
		}
		else if (leftActiveEdge->leftArc->GetP(d) > FLT_EPSILON)
		{
			resPoints[0].push_back(leftActiveEdge->leftArc->GetPointAtAngleParam(d, this->GetRotateAngle()));		
		}
		else if (rightActiveEdge->rightArc->GetP(d) > FLT_EPSILON)
		{
			resPoints[0].push_back(rightActiveEdge->rightArc->GetPointAtAngleParam(d, this->GetRotateAngle()));		
		}
		else
		{
			resPoints[0].push_back(majorFocus);
		}
		return resPoints;
	}
	float2 angles = GetAngles(d);
	float2 angles2 = make_float2(1, -1);
	if (angles.x < -10 || angles.y > 10)
	{
		assert(false);
	}
	if (maxAngle <= M_PI && angles.y <= M_PI)
	{
		angles.y = min(angles.y, maxAngle);
		angles.x = max(angles.x, minAngle);
	}
	else if (maxAngle > M_PI && angles.y <= M_PI)
	{
		angles2 = make_float2(angles.x, maxAngle - M_PI);
		angles.x = minAngle;
	}
	else if (maxAngle <= M_PI && angles.y > M_PI)
	{
		angles2 = make_float2(angles.x, maxAngle);
		angles.x = minAngle;
		angles.y -= 2 * M_PI;
	}
	else
	{
		angles.x = max(angles.x, minAngle);
		angles.y = min(angles.y, maxAngle);
	}

	float deta = (maxAngle - minAngle) / 360;
	if (angles.x < angles.y)
	{
		resPoints.push_back(std::vector<float2>());
		int pos = resPoints.size() - 1;
		for (float angle = angles.x; angle <= angles.y; angle += deta)
		{
			resPoints[pos].push_back(GetPointAtAngleParam(d, angle));
		}
		resPoints[pos].push_back(GetPointAtAngleParam(d, angles.y));
	}


	if (angles2.x < angles2.y)
	{
		resPoints.push_back(std::vector<float2>());
		int pos = resPoints.size() - 1;
		for (float angle = angles2.x; angle <= angles2.y; angle += deta)
		{
			resPoints[pos].push_back(GetPointAtAngleParam(d, angle));
		}
		resPoints[pos].push_back(GetPointAtAngleParam(d, angles2.y));
	}
	/*TRACE("!!!!!!!!!!!!!!!!!!!!!!!!!\n");
	TRACE("resPoints  size %d ", resPoints.size() );
	for(int i = 0; i < resPoints.size(); ++i){
	TRACE("resPoints %d is %d\n: " , i , resPoints[i].size());
	}
	TRACE("\n");*/
	return resPoints;
}







struct __align__(8) ArcHeapEle{
	float angle;
	EllipseArc* arc;
	__device__ ArcHeapEle(float _angle,EllipseArc* _arc){
		angle = _angle;
		arc = _arc;
	}
	__device__ bool operator<(const ArcHeapEle& other) const{
		return angle < other.angle;
	}
};

typedef EllipseArc* pEllipseArc;

template<class T>
struct myPoint{
private:
	T* h_point;
	T* d_point;
	int* h_length;
	int* d_length;
	int sizePerThread;
	int totalThreadNum;
	int totalSize;
public:
	T* hostPoint(){
		return h_point;
	}
	int* hostLength(){
		return h_length;
	}
	T* devicePoint(){
		return d_point;
	}
	int* deviceLength(){
		return d_length;
	}
	void initial(const int& _sizePerThread,const int& _totalThreadNum){
		sizePerThread = _sizePerThread;
		totalThreadNum = _totalThreadNum;
		totalSize = sizePerThread * totalThreadNum * sizeof(T);
		allocHost();
		allocDevice();
		for(int i = 0; i < totalThreadNum; ++i){
			h_length[i] = 0;
		}
		for(int i = 0; i < sizePerThread*totalThreadNum;++i){
			//h_point[i].initial();
		}
	}

	void insert(const int& threadPos,const T& tmpElement){
		h_point[threadPos * sizePerThread + h_length[threadPos] ] = tmpElement;
		h_length[threadPos]++;
	}

	void memcpyH2D(){
		cutilSafeCall( cudaMemcpy( d_point,h_point,totalSize,cudaMemcpyHostToDevice));
		cutilSafeCall( cudaMemcpy( d_length,h_length,sizeof(int)*totalThreadNum,cudaMemcpyHostToDevice));
	}

	void memcpyD2H(){
		cutilSafeCall( cudaMemcpy(h_point,d_point,totalSize,cudaMemcpyDeviceToHost));
		cutilSafeCall( cudaMemcpy(h_length,d_length,sizeof(int)*totalThreadNum,cudaMemcpyDeviceToHost));
	}
	T& atHost(const int& threadPos, const int& pos){
		return h_point[threadPos*sizePerThread + pos];
	}
	int& length(const int& threadPos){
		return h_length[threadPos]; 
	}
	void unalloc(){
		if( h_point != NULL ){
			free(h_point);
		}
		if( h_length != NULL ){
			free( h_length );
		}
		if( d_point != NULL )
			cudaFree( d_point );
		if( d_length != NULL )
			cudaFree( d_length );
	}
private:
	void allocHost(){
		//h_point = (T*)malloc(sizePerThread*totalThreadNum*sizeof(T));
		//h_length=(int*)malloc(sizeof(int)*totalThreadNum);
		cutilSafeCall( cudaHostAlloc((void**)&h_point,sizePerThread*totalThreadNum*sizeof(T),cudaHostAllocDefault));
		cutilSafeCall( cudaHostAlloc((void**)&h_length,sizeof(int)*totalThreadNum,cudaHostAllocDefault));

	}

	void allocDevice(){
		cutilSafeCall( cudaMalloc((void**)&d_point,totalSize));
		cutilSafeCall( cudaMalloc((void**)&d_length,sizeof(int)*totalThreadNum));
	}

};


//float2* h_sites;
//float2* d_sites;
//int* h_sitesLength;
//int* d_sitesLength;
myPoint<float2>sites;

//int* h_vertexPoolLength;
//float2* h_vertexPool;
//int* d_vertexPoolLength;
//float2* d_vertexPool;
myPoint<float2>vertices;

//int* h_edgePoolLength;
//int2* h_edgePool;
//int* d_edgePoolLength;
//int2* d_edgePool;
myPoint<int2>edges;

myPoint<EllipseArc>arcs;
myPoint<EllipseArcNoPt>arcsNoPt;

myPoint<float2>bounds;

myPoint<int>iteratorItmes;

//ActiveEdge* d_activeEdgePool;
//int* d_activeEdgePoolLength;
myPoint<ActiveEdge>activeEdges;
myPoint<ActiveEdgeNoPt>activeEdgesNoPt;



float radiusVaring;
//float2* h_center;
//float2* d_center;
myPoint<float2>centers;
myPoint<float2>distances;
myPoint<int>headArcsPos;

myPoint<int>iterateTimes;



int threadNum;
int blockNum;

myPoint<Event> eventHeaps;
myPoint<EventNoPt> eventHeapsNoPt;
myPoint<ArcHeapEle> arcHeaps;


float tmpDisFromCenterScale;

#endif

